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Abstract:  Formulation  is  given  for  efficient  parabolic  equation  solution  of  radiowave 
propagation  in  inhomogeneous  atmosphere  and  over  irregular  terrain.  Both  standard 
and  wide  angle  parabolic  equation  derivations  are  presented.  Impedance  boundary 
conditions  are  used  to  characterize  the  ground.  A  tropospheric  boundary  condition 
based  on  the  exact  solution  of  Schrodinger  equation  in  a  quarter  plane  is  derived.  To 
permit  efficient  modeling  of  the  irregular  boundary,  the  parabolic  equation  together 
with  the  boundary  conditions  is  transformed  into  a  numerically  generated  curvilinear 
coordinate  system.  Finally,  formulation  is  presented  for  a  finite  difference  solution 
using  Crank- Nicolson  implicit  scheme. 


1.      Introduction 

It  is  well  known  that  multipath  fading  can  significantly  effect  the  link  reliability  in  a 
communications  system  or  target  detectability  in  the  case  of  radar.  The  path  between 
a  transmitted  and  a  receiver  is  often  obstructed  by  natrual  or  man-made  obstacles 
such  as  hills,  buildings,  atmospheric  layers,  trees,  rain,  fog,  etc.  Propagation  outage 
due  to  multipath  fading  depends  in  a  complicated  manner  on  propagation  climate, 
terrain  features,  path  length,  radio  frequency,  and  fade  margin.  In  the  case  of  atmo- 
spheric multipath  fading,  interference  due  to  two  or  more  super- refracted  rays  arriving 
at  the  receiver  via  different  paths  can  lead  to  a  complete  loss  of  signal.  Other  patho- 
logical phenomenon  such  as  obstruction  fading  (caused  by  sub-refractive  atmospheric 
effects)  and  ducting  (caused  by  extreme  super- refractive  effects)  are  also  possible. 
Such  phenomenon  are  more  common  in  warm  and  tropical  climates,  particularly  near 
shores,  where  elevated  inversions  are  formed  easily  due  to  the  large  temperature  and 
partial  pressure  differentials.  Reflection  multipath  fading,  which  is  due  to  interfer- 
ence between  the  direct  and  the  ground  reflected  ray  depends  strongly  on  the  terrain 
geometry  and  ground  constants.  Moreover,  elevated  terrain  features  could  completely 
mask  a  receiver  from  a  transmitter  leading  to  severe  loss  of  signal  (in  some  cases,  it 
is  advantageous  to  site  antennas  behind  hills  to  provide  shielding  against  undesirable 
interference).  It  is  very  important  to  assess  the  effects  of  environment  on  the  link.  A 
computer  model  that  can  take  into  account  a  given  refractive  index  profile,  terrain 
elevation  data,  and  varying  ground  parameters  will  be  very  helpful  in  predicting  the 
link  performance. 

In  this  report  we  present  formulation  details  for  an  efficient  numerical  solution  of 
wave  propagation  in  an  inhomogeneous  atmosphere  and  over  irregular  terrain  using 
parabolic  equation. 

Unlike  all  other  previous  formulations  of  the  parabolic  equation,  we  will  use  a  modified 
Helmholz  equation  for  propagation  in  an  inhomogeneous  atmosphere  as  suggested  by 
Maxwell's  equations  (all  previous  formulations  use  a  Helmholtz  equation  which  is 
only  true  for  fields  in  a  homogeneous  medium).  Because  the  parabolic  equation  is  a 
full-wave  method,  it  will  include  all  aspects  of  wave  propagation  such  as  reflection, 
refraction,  diffraction,  and  surface  wave  propagation.  In  this  respect  it  is  far  superior 
to  the  commonly  used  ray  method. 

Parabolic  equation  approximation  to  an  elliptic  partial  differential  equation,  which 
the  true  fields  satisfy,  has  proven  to  be  a  viable  approach  for  studying  propagation 
problems  in  underwater  acoustics.  The  method  is  just  gaining  popularity  with  the 
electromagnetic  community.  Although  the  parabolic  equation  regards  waves  as  es- 
sentially traveling  one-way,  it  allows  a  rapid  solution  of  the  fields  by  way  of  marching 
along  the  range  starting  from  an  initial  range.  Another  advantage  of  the  PE  method 


compared  to  the  ray  methods  is  that  it  is  valid  even  in  the  shadow  region  where  the 
simple  ray  methods  completely  break  down.  Furthermore,  it  appears  to  be  the  only 
practical  method  for  predicting  propagation  over  long  ranges  (greater  than  1  km)  over 
a  wideband  from  HF  (a  few  MHz)  through  SHF  (a  few  tens  of  GHz).  The  method  is, 
however,  not  without  limitations.  In  its  standard  form,  the  accuracy  of  the  method 
is  limited  to  waves  traveling  essentially  within  ±10°  from  horizontal.  Furthermore, 
treatment  of  the  boundary  conditions  on  the  uneven  terrain  is  difficult. 

The  method  we  propose  to  use  will  attempt  to  remove  both  of  these  deficiencies. 
Firstly,  we  use  a  Helmholtz-like  elliptic  equation  to  describe  the  fields  and  arrive 
at  a  wide-angle  parabolic  equation  subject  to  certain  approximations.  To  facilitate 
imposition  of  the  boundary  conditions  on  the  irregular  terrain,  the  equations  will  be 
transformed  to  a  body-fitted  curvilinear  coordinate  system.  The  PE  will  be  solved 
using  finite  differences  on  a  non-rectangular  mesh.  This  is  a  major  departure  from 
previous  approaches  which  will  not  only  make  the  method  more  efficient  but  also 
more  accurate. 

In  Section  2,  we  derive  the  exact  equations  satisfied  by  2D  fields  in  an  inhomogeneous 
atmosphere.  Impedance  boundary  conditions  are  used  to  characterize  the  ground.  In 
Section  3,  we  present  details  on  the  impedance  boundary  conditions.  Starting  from 
the  exact  equations  presented  in  Section  2  for  the  fields,  we  derive,  in  Section  4,  a 
parabolic  equation  (PE)  valid  for  narrow  angle  propagation.  This  case  is  termed  as 
the  standard  PE.  The  standard  PE  is  generally  valid  for  propagation  angles  that 
are  within  ±10°  from  horizontal.  To  accomodate  waves  traveling  at  larger  angles,  we 
present  the  derivation  of  a  wide  angle  PE  in  Section  5.  To  truncate  the  computational 
domain  we  derive  boundary  conditions  on  an  upper  boundary,  which  are  termed  as 
the  tropospheric  boundary  conditions.  The  derivation  is  based  on  the  solution  of 
Schrodinger  type  parabolic  equation  in  a  quarter  plane  x  >  0,  y  >  0.  This  is  presented 
in  Section  6. 

For  an  efficient  numerical  implementation,  we  transform  the  differential  equation  and 
boundary  conditions  to  a  curvilinear  coordinate  system.  This  is  presented  in  Section 
7.  Details  on  the  numerical  generation  of  the  curvilinear  coordinate  system  are  given 
in  Section  8.  Finally,  in  Section  9,  we  present  steps  leading  to  a  finite  difference 
solution  of  the  equations  using  a  Crank-Nicholson  implicit  scheme. 


2.      Solution  of  2D  Fields  in  an  Inhomogeneous  Medium 

Consider  an  electric  source  producing  fields  in  an  inhomogeneous  region  as  shown 
in  the  figure  below.  Let  us  assume  that  both  the  sources  and  the  medium  are  two- 
dimensional  in  nature  in  that  all  quantities  are  independent  of  the  ^-coordinate.  As 
in  the  case  of  a  homogeneous  medium  the  fields  can  be  decomposed  into  a  TEZ  case 
(vertical  polarization)  and  a  TMZ  case  (horizontal  polarization).  It  is  assumed  that 
propagation  takes  place  in  the  x?/-plane. 
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From  Maxwell's  equations,  we  have 

V  x  E    =    —jujfiH 

V  x  H    =    jujtE  +  J 


(1) 
(2) 


In  the  case  of  vertical  polarization,  the  fields  could  be  written  in  terms  of  the  z- 
component  of  the  magnetic  field,  H  =  zHz,  {TEZ  fields).  Substituting  into  (2)  we 
have 


V  x  H  =  V  x  {zHz)  =  S7HZ  x  z  =  jueE  +  J  =»  tE 
In  a  source-free  environment,  we  have 


VHZ  x  z  -  J 


cE  = 


VHZ  x  z 


Substituting  into  (1)  we  get 


1  „      fVHz        \           z  „    fVHt\ 
V  x  E  =  —V  x x  z )  =  -—V  • =    -jufiH 

JU  V       t  J  JU!  \       € 


The  equation  satisfied  by  the  magnetic  field  is  then 

iu        \     e    ) 


v.p^)+wvr,  =  o 


Vertical  Polarization 


Once  the  magnetic  field  is  determined  the  electric  field  is  given  by 


E  = 


zxVHz 

jut 


(3) 


(4) 


Note  that  Hz  does  not  satisfy  the  Helmholtz  equation  unless  t  is  constant. 
For  horizontal  polarization  on  a  similar  analysis  with  E  —  zEz  shows  that 

Horizontal  Polarization 


(5) 


If  the  medium  is  non-magnetic,  //  =  fi0  and  Ez  satisfies  the  Helmholtz  equation.  We 
will  characterize  the  ground  in  terms  of  impedance  boundary  conditions  [3]. 

We  may  combine  the  vertical  and  horizontal  cases  shown  in  (3)  and  (5)  into  an 
equation  of  the  form 

V  •  (aVVO  +  H  =  0  (6) 


where 


a    =     < 


-        TE  or  Vertical  Pol. 
e 

TM  or  Horizontal  Pol. 


(3    =    akln2(x,y)  >  0 
Hz       TE  Pol. 


0    = 


(7) 

(8) 
(9) 


E2 


TM  Pol. 


The  quantity  n(x,y)  is  a  position  dependent  refractive  index  of  the  medium.    The 
partial  differential  equation  (PDE)  given  in  (6)  is  elliptic,  for  we  have 

d_(    d±\        d_  (    df\ 
dx  \    dx  J       dy  \    dy  I 


or 


Q 


<920      d2iP 
+ 


dx 


dr 


dib  da      dib  da 
ox   ox       oy    oy 


The  discriminant  for  the  above  PDE  is  —  1  <  0  implying  elliptic  nature  of  the  equation. 
Equation  (6)  could  be  expressed  as 


W  +  -    7T V'x  +  "    ^y  +  ~ll>  =  0 

a   ox  a   oy  a 


For  a  non-magnetic,  loss-less  medium,  we  have 

1 


) 

Vertical  Pol. 

c0er 

a  =  < 

1 

• 

Horizontal  Pol. 

.     //o 

polarization 

1    da                5/1 

:) 

1    der            1    dn2 

a   dx 

r  dx  \c, 

tT    dx           n2    dx 

Similarly, 
Letting 


2  dn     ,    .     ,        r       .      .    ,     , 

= — —     [n  is  the  refractive  index) 

n  dx 


I    da  2  dn 

a  dy  n  dy 


"i{x,y)    =    < 


a2{x,y)     =     < 


2   dn 
n   dx 

0 

2  dn 
n  dy 
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Vertical  Pol. 

Horizontal  Pol. 
Vertical  Pol. 

Horizontal  Pol. 


equation  (10)  may  be  expressed  in  the  form 


(10) 


(11) 


3.      Impedance  Boundary  Condition 

Impedance  boundary  condition  relates  the  tangential  components  of  electric  and 
magnetic  fields  at  the  interface  of  two  media.  If  fr  is  a  unit  normal  and  s  is  a  unit 
tangent  as  shown  in  the  figure  below,  the  boundary  conditions  are  given  by  [3] 


s 


£ 


\ 


V 


s  =  x  cos  9  +  y  sin  9 
v  =  —x  sin  9  -f  y  cos  9 

x  =  s  cos  0  —  if  sin  0 
y  =  5  cos  #  +  if  cos  0 


h(h£)=  -r/0As/>  x  #  (12) 

where  As  =  Zs/t)0  is  the  surface  impedance  normalized  to  the  free  space  value  t/0. 
The  equation  may  also  be  expressed  as 

[y  ■  E)v  —  E  =  —t)0Asu  x  //  =>  0  x  E    =    t}0AsC>  x  (/>  x  //) 

-    itoA.  [(*  ■  J?)*  -  if]  (13) 

The  surface  impedance  is  determined  from  the  intrinsic  impedance  of  the  medium 
by  considering  plane  wave  reflections  from  the  interface.  The  complex  propagation 
constants,  7l5  72,  and  the  intrinsic  impedances  ^  and  r/2  in  terms  of  the  media 
constants  are  indicated  in  the  figure. 


7{      =      jUflrifi0((Ti  +  jWCoCrl) 


—     —k0^rl  I  erl  —  j 


€v<€o,/xr(//0 ,  <T( 


.  o-i 


we0 


T—> 7 >        >        /        f 


— =  Vo\  


V2 


Voy 


1Ht2 
Cr2 


Vo  = 


According  to  Snell's  law,  we  have  7j  sin  9X  =  ~f2  sm  @t 


The  plane  wave  reflection  coefficients  for  the  vertical  and  horizontal  polarizations,  Rv 
and  Rh  are  [7] 

772  sec  9t  —  ■q1  sec  6t 


Rh 
Rv 


Tj2  sec  0t  +  77j  sec  62 

T}2  COS  6t  —  T]i  cos  9t 
Tj2  COS  6t  —  T]i  COS  0t 


Surface  Impedance  Zs    —  r/2  sec  6t 

ZVS    =  T}2  cos  9t 


ah  V2  sec  6t  __  1J2 


V\ 
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1  -  (  —  sin  Oi 


and 


aY  = 


^rl^rcl  2    / 

-  cos  xpi 


Vt2(-tc2 


1-1/2 


1/2 


^t2^tc\ 


V   I1t\^tc2 

For  the  special  case  of  //1  =  /z2  =  A*o,  Ci  =  0, 


1 COS     Ipi 

Ht2(-tc2 


-.1/2 


A?    = 
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Cr2  -JOV2 
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£rl 


1-1/2 


Crl 


V  Cr2   -  JCTr2 

For  normal  incidence  ipt  =  90° 

A"  =  Av  = 

s  s 


Cr-2  -jCTr2 

(■rl 
f-r2  -  ]<?t2 

Crl 


cos2  Ipi 
cos2  ^,- 


11/2 


(14) 
(15) 


V   tr2  -  J°t2 

For  the  2-D  case,  impedance  boundary  conditions  for  vertical  polarization  can  be 
simplified  as 

Ox  E  =  t]0A^  [{0  ■  H)0  -H]=  -Vo^H 

Taking  a  dot  product  on  both  sides  with  z 
z-(0xE)  =  -Vo^Hz 

Since  z  x  0  —  —  s,  we  have 
-s-E  =  -Vo^Hz 

Substituting  from  equation  (4)  for  E,  we  get 

zxVH 


JUJC 


=  -voa:h 


s±J-z 


i>VH,  =  juctioA;Hz  =  jk0tTAvsH 


S    "Z 


i.e.,  the  above  can  be  put  conveniently  in  the  form 

IBC  Vertical  Pol. 


dHx 


dv 


-jkoerA^Hz  =  0 


A  similar  analysis  for  horizontal  polarization  yields 


dE'       -u        l    r       n 


IBC  Horizontal  Pol. 


This  may  also  be  obtained  by  resorting  to  duality. 
For  a  perfectly  conducting  material,  we  have  A^v  =  0  and 

dH7 


dv 


=    0         Vertical  Pol. 


E7     —    0         Horizontal  Pol. 


(16) 


(17) 


4.      Standard  PE  Derivation 


We  will  make  some  approximations  and  cast  (11)  in  the  form  of  a  parabolic  equa- 
tion which  permits  a  rapid  numerical  solution. 


Let  4>{x,y) 


,-jko* 


y/x 


■u(x,y) 


Then 


*    =     (+ -*»-&) 


-jk0x 


2x3/2J        y/i~ 


(«i  -  jk0u)- 


,-Jk0x 


y/x 


X  — >  oo 


,-jk0x 


*Py     —     uy — 7=~  ■>      Uyx  =  uTy  ~  (uxy  —  jk0uy)- 


y/x 


,-jkox 


y/x 


^yy     =     ui 


,-jkox 


Substituting  into  (17)  we  get 


,      uxx  ~  (uxx  -  2jk0ux  -  k0u)- 


,-jk0x 


y/x 


uxx  +  uyy  —  2jk0ux  -f  a\Ux  —  2jkQa-[u  -f  a2uy  -\-  (n2  —  1)A:qW  =  0 

or 

Wu  +  uyy  +  (ai  -  2jk0)ux  +  a2uy  +  (n2  -  1  -  2j  —  J  A;Ju  =  0 

If  now  we  impose  the  approximation  that 

1  /9 

\uxx\  <  (a J  +  4ArJJ       \ux\  , 


we  obtain 


1 


Ur    = 


(aj  -  2;A:0) 


uyy  +  a2uy  +  (rc2  -  1  -  2j  —  j  Arjj  u 


or 


(18) 


This  is  the  exact  form  of  narrow  angle  PE  approximation.    We  would  also  like  to 
express  the  impedance  boundary  condition  in  terms  of  the  'u'  functions. 


S 


x    =  s  cos  9  —  v  sin  6 

dx 

xu  =  — —     =  v  ■  Vx  =  v  •  x 
ov 

=  —  sin  6 

v    =  —x  sin  6  +  y  cos  6 


10 


For  vertical  polarization  we  have 
dH 


di 


-jkoerAv$Hz  =  0 


dHz 
du 


-jkon2A^Hz  =  0 


Hg(x,y)    = 


y/x 


■u(x,y) 


dHz 

du 


—       I  Uv  —  JKqX^U  — 


ux, 


,—jkox 


2x3/2J      yfi 


{uv  —  jk0x„u) 


,-jkox 


noo 


u>„  -  jk0(n2Avs  +  x„)u  =  0 


y/x 

IBC  Vert.  Pol. 


Similarly  for  horizontal  polarization,  we  have 

IBC  Horz.  Pol. 

n  '    i 

s  / 

We  combine  the  two  by  defining 

'  -jk0{n2Avs  -sin0)        Vert. 
c\  =  < 


~jk0  (  -£jj  -  s'mOj        Horz. 


and  writing  as 


uu  +  c-[U  =  0 


(19) 


The  parabolic  equation  given  in  (18)  is  valid  for  propagating  angles  close  to  horizon- 
tal (±10°  in  practice)  [1].  To  accomodate  waves  at  higher  angles  we  would  need  a 
wide  angle  parabolic  equation  whose  derivation  is  accomplished  through  a  pseudo- 
differential  operator  formalism  [1]. 
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5.     Wide  Angle  PE  Derivation 


Let  us  assume  that  the  media  constants  are  independent  of  horizontal  range  so  that 
a(x,y)  =  a(y),  n(x,y)  =  n(y).  We  then  have  from  (17) 

urx  +  uyy  -  j2k0ux  4-  -  — uy  +  {n2{y)  -  l)fcju  =  0 

a    ay 


Let 


P 


d_ 
dx 


Q  = 


l  d2 


+ 


1    da  d 


\  kl  dy2       kla  dy  dy 


+  n' 


Pseudo  differential 
operator  in  y 


(20) 


Using  this  notation,  the  above  PDE  may  be  expressed  as 

P2-2jkoP  +  (Q2-\)k2o\u  =  0 

which  we  may  factorize  as 

(P  ~  jk0  -  jk0Q)(P  -  jk0  +  jkoQ)u  =  0 

To  see  this  we  expand  the  operator  on  the  left  hand  side  to  get 

P2  ~  jk0P  -  jkoQP  -  jkQP  -  k20  -  k20Q 
+  jk0PQ  +  k20Q2  +  k2Q 

Since  PQ  =  QP,  we  have  the  desired  result.  The  first  operator  in  (20)  denotes 
an  incoming  wave  (w.r.t.  x)  and  the  second  an  outgoing  wave.  We  retain  only  the 
outgoing  wave  to  obtain 

Pu=  -jk0(Q-\)u  (21) 

The  square  root  operator  Q  is  global  in  nature  and  we  would  like  to  make  some 
approximations  to  derive  a  local  operator  from  it.  (It  is  global  because  when  expanded 
in  terms  of  series,  it  will  contain  terms  of  all  derivatives).  Now 


Q  = 


kl  a   dy   dy       kl   dy2 


1/2 


Let  us  rewrite  Q  as 


Q  = 


l  +  [n2(y)-l]  +  - 


1       da         d 


+ 


cP 


<C1   normally 


a  d(k0y)d(k0y)       d{k0y)2 


iall 


1/2 
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which  may  be  expressed  as 
where 


q  =  vTTV 


2      n         1     da    d        Id2 

V    =    n    —  1  H 1 

kla  dy   dy       k%  dy2 

=    a  small  operator! 

Treating  V  as  an  algebraic  factor  <  1,  we  may  derive  the  following  rational  approxi- 
mation (pade(l,  1)) 

—       i  +  fv 

j —     (Claerbout) 

4  +  3V  2V 

4  +  V        1  +  4  + V 

2V 


g  =  vTTV 


so  that 


Q-l  = 


4  +  V 


Substituting  into  (21)  we  arrive  at 

Pu  =  jk0{Q-  l)u=>  Pu 


-2jk0V 

4  +  V 


or 


(4  +  V)Pu  =  -2jk0Vu 


tt  •  r  I  da 

Using  the  fact  that  a2  =  —  -7—,  we  get 

a   dy 


,  ,      MI,  5        d2 


«i  =  -2jk0 


d        d: 


{nz  -  l)k'0  +  a2—  + 


dy     cty2 


(22) 


The  above  equation  is  a  wide-angle  parabolic  equation  valid  for  propagation  up  to 
±20°  [1].  Other  approximations  could  be  obtained  by  considering  higher  order  pade 
approximants. 
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6.      Boundary  Condition  on  the  Upper  Boundary 

To  truncate  the  computational  domain,  we  consider  a  point  high  enough  where  the 
atmosphere  is  homogeneous  with  n  —  1.  The  governing  equation  in  the  homogeneous 
region  becomes 

J 
Ux  ~     2k0Uyy 

Let  us  derive  boundary  conditions  on  a  horizontal  interface  y  =  y0.  For  the  sake 
of  simplicity  we  will  derive  boundary  conditions  of  the  mixed  type  on  the  interface 
y  =  0  (instead  of  y  =  t/o)  given  the  initial  data  on  the  line  x  —  0  and  boundary  data 
on  the  line  y  =  0.  Our  derivation  is  based  on  the  use  of  Fourier  sine  transforms  as 
suggested  in  [2].  Although  the  basic  philosophy  of  our  approach  coincides  with  that 
in  [4],  some  of  the  details  and  the  final  results  are  slightly  different  from  the  latter. 

Consider  the  parabolic  equation  uyy  —  2jkux  =  0  in  a  homogeneous  region  x  >  0, 
y  >  0,  where  A;  is  a  complex  constant, 


A* 

;^u=f(v) 

/ 

Uyy      " 

2 

•  ix-  ?(«) 

-       / 

> 7- I 7 7 

"7 7 

— 7— 

subject  to  the  initial  condition  u(0,  y)  =  /(y),  0  <  y  <  oo,  and  the  boundary  condition 
u(x,0)  =  g(x),  0  <  x  <  oo.  We  assume  that  u(x,oo)  — >  0  and  uy(x,oc)  — >  0.  The 
equation 


Lyy 


=  2jkui 


(23) 


is  of  Schrodinger's  type.  We  will  treat  the  lossless  case  having  a  real  value  of  k  as 
the  limiting  case  of  the  lossy  problem  having  k  =  k0  —  je,  t  >  0.  We  will  solve  the 
problem  using  Fourier  sine  integral. 


Let 


Us{x,  X)  =  J-  I     u(x,  y)  sin  Ay  dy 


(24) 


\\ 


Then  using  integration  by  parts,  we  see  that 

-/     uyy(x,y)s'm\ydy    =     J-  I  uy(x,  y)  sin  Xy 


y=0 


y=o 


/•oo 

A    /      u(x,  y)  sin  Xy  dy 


—Xu(x,  y)  cos  Xy 
Because  uy(x,  oo)  — ►  0  and  u(x,  oo)  — ►  0,  we  have 

y~  /     uyy(x,y)  sm  Xudy    =     J-Xu(x,  0)  -  X2Us{x,  X) 


■Xg(x)-X2Us(x,X) 


(25) 


Multiplying  both  sides  of  (23)  with  w2/7rsinAy,  integrating  over  y  =  0  to  oo  and 
making  use  of  (25),  we  have 


or 


'-\g(x)  -  X2Us(x,  X)  =  2jk^-Ua(x,  X) 

IT  OX 


lu,(x,\)  +  £-kU,(z,\)  =  ±Jl9(x) 


We  may  rewrite  the  above  equation  as 
d 


d: 


X       2 


Us(x,X)e^/2^    =  JL.J±g(x)elx2Wk)* 


(26) 


2  j  k  V  7T  " 

Replacing  the  dummy  variable  x  in  (26)  with  r  and  integrating  both  sides  over  r  =  0 
to  x,  we  arrive  at 


Us{r,X)e^'2^T 


x     2  r 


-  I'    g(r)e^l2^Tdr 

TT  Jt=0 


or 


But 


r=0  W 

U.{zt  X)  =  U.(0,  \)e-l*W*  +  -^  f   g(T)e*2IW-x)dr 

Zj  K  V   7T  7t=0 

/2   y00 

t/s(0,A)     =     W-yo     u(0,y)sinAj/c/y 

=  v~ J0  f(y) sin A^ dy =  F*(x) 
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Ua(x,  A)  =  Fa(A)e-^2/2^  +  £r\[lr   9(t)** ***-*&  (27) 

ZjK  V   7T  Jt=0 

Finally,  taking  the  inverse  sine  transform  on  (27)  we  get 
u(x,y)    =     \   -  /     sin  Ay  Us(x,  \)d\ 

V  7T  JO 

=     iFf  Fs(A)e-(A2/2^sinAyrfA 

V  7T  J\=0 

+  -^7  /°°  A  sin  Ay  f    y(r)e(A2/2^T-rW  rfA 

7T  ZJK  J\=0  Jr=0 

0  roo        roo 

=     -/        /      /(r)sinArrfre-(A  /2^)xsinAydA 

7T  7A=0  ./t=0 

+--V  f    <7(r)  /°°  AsinAye^^-^A^r 
7T  2jK  Jt=0  Jx-o 

=     -  /       f(T)  /       [cos(T  -  ?/)A  "  cos(r  +  V) Al  e_(    /2J  )xdA  dr 
it  JT=o         jk  \    oy  J  J\=o 

1  f<X>  TOO  „ 

=     -/      f(T)  [cos(y-r)A-cos(y  +  r)A]e-(A  /2]k)xdXdr 

IT  Jt=0  J\=0 

~  f    y(r)|-{  rcosyAe-(A2/2^-^A}rfr 
7rj A:  7t=o  oy  (.Jo  ) 

Defining 

roo 

A'(x,y;  x0,y0)  =  /     cos(y  -  y0)Ae_(A  /2^)(l"Xo)(fA,      for  x  >  x0, 
Jo 

we  write  the  expression  for  u(x,y)  as 

1    /"°° 
u(x,y)     =     -  /(r)[A'(i,y;  0,r)-  A'(x,y;  0,-t)]<*t 

7T  7t=0 

-^jC.*^™^*  (28) 

We  now  evaluate  the  integral  for  K .  Consider 

1(a)     =     [°°  cosa\e-{x2/2jk){T-xo)d\,  x  >  x0 

J\=o 
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Jx=o 

roo 

Jx=o 


Because  of  the  exponential  decay,  we  may  differentiate  under  the  integral  sign  to 
obtain 


da 


1(a)    =     -  f°°  XsmaXe-^^W^-^dX 
Jx=o 


=      /       sinaA — -t- 
J\=o  [x 


■(A2(x-xo)/2|/c|2)(£-jA:o) 


(x  -  x0)(e  -  jk0)/\k\2 
sinaAe-A2(l-Xo)(£-j/co)/2|fc|2 


dX 


(x  -  x0)(e  -  jk0)/\k\2 


a 


k? 


^=0      {x  -  x0){e  -  jk0) 


roo 

/      cos  a 
Jo 


Ae-(A2(x-x0)(<-jfco))/2|M2dA 


ja\kf 


{X  —  XqJK 


roo 

-  /      cos  ay 
'  Jo 


e-X>(x-I0)/2jkdX  =  _ 


jak 


(x  -  x0) 


1(a) 


dl(a)  jak 

—  /(a)  =  0 


d_ 

da 


da  (x  —  Xo) 

Now  /(a  =  0)  can  be  written  as 

l(a  =  0)=   re-(A2(x-,o)/2N2)(^o)^A 

Jo 

Let  us  view  this  integral  in  the  complex  A-plane. 

I~0H         -     —       "       e 


rjrv 


H(A>>  -    T       2*, 


4 


Cffvv^-^^"- 


AvgU)  =  f| 


v  R*(A) 


Complex  A-Plane 


(29) 
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For  the  integral  to  converge  for  (x  —  x0)  >  0,  Re[X2(e  —  jk0)}  >  0,  i.e., 

Re[X2(-jk')}  >  0  =»  Im{X2km)  >  0 


or 


or 


or 


0  <  Arg(A2F)  <  7T 


0  <  2Arg(A)-  Arg(fc)  <  tt 


Now 


1  7T  1 

-  Arg  {k)  <  Arg  ( A)  <  -  +  -  Arg  (k) 
Arg(A;)  =  -tan-1(p)«-^,  i- <  1 

\  /C0  /  /Cq  /Cq 

_1_       A      (X)       -      — 
Zkq  Z        Z,Kq 

I(a  =  Q)  =  J  e-x2(x-x°WVW2d\ 

where  C  lies  in  the  region  of  convergence. 

In  particular,  let  us  choose  a  line  from  0  to  oo  along  the  line  Cw  defined  by 


A  = 


'   2jk 

X  —  Xq 


■w  w  —  0  to  oo 


I-vw(A)  I 


Av^(A)  - 


TV 


(30) 


(31 


>  Re  (A) 


On  this  path, 


I{a  =  0)  = 


Arg(A)  =  I  +  -ArgW  =  I-- 


2jk 


{X  —  Xq)  J w 


/•oo 

Jw=0 


'         7TJk 

2(x  -  x0) 
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1(a)  = 


TTjk 


2(x  -  x0) 


,-ja2k/(2(x-x0)) 


K(x,y;  x0,y0)  =  J     *2*      e-««.JVW«. » 

V  2(x  -  xq) 


(32) 


Substituting  this  into  (28),  the  field  u(x,y)  for  the  initial-boundary  value  problem  is 
given  by 

7T  Jt=0  V     2x      *•  J 


1 

TTJk 


It  is  easy  to  see  that 

dK 


£    =    /;-^costo-?0)A]e-^-^»,A 


2f/ 


52/\ 


5r 


too 

/     -A2cos[(y-y0)A]e-A^-T^2^<fA 


and,  so  for  x  ^  x0 


d2K     2-kdK  dK_        1    d2K 

dy2  dx  dx        2jk  dy2 


which  is  the  same  equation  satisfied  by  u.  Now 

u(x,y)     =     JlL  f°°  f(T)\e-^- 
V  lirx  Jt-o  l 


02/(2x)  _  e-jk(y  +  ry/(2x) 


1      fx  d 

T  /      g(T)—K{x,y;  r,0)c?r 

7tjk  Jt=o         oy 


du 
oy 


y->0+ 


2irx  Jo  x 

vrkLg{T)$K{x'y''T'0) 

— ;  f  g(T)2jk—K(x,y-  r,0) 

IT  IK  J0  OX 


dr 


dr 


y-*0+ 


(33) 


(34) 
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in  view  of  (34).  Furthermore,  we  note  that 

—K(x,y;x0,y0)  =  -  —  K(x0y;x0yo) 
ox  oxq 

Using  this  and  evaluating  the  first  integral  by  parts,  we  get 


9     (        \ 

oy 


♦0+ 


V    TTX      {  t=q         Jo 

2    [x  d 

+  "  /      g{T)-^-K(x,y\  r,0)  dr 

TV  Jr=0  OT  n. 


f2jk 


—  f-/(0)-  r  fr(r)e-^^dr} 

TTX     I  JO  J 


2 

+- 

7T 


g(T)k(x,y,T,0) 


T_0    "  f    gT(r)K(x,0-  r,0)di 

„  +  Jr=0 

y  —  Cm 


=     +V-^°)  +  V-/    Mr)*-**™* 

V     TTX  \     TTX    Jo 

Jh. 


2    <c\\  rjk  2lk   [X      9r{r), 


2x 


=     J— If  (0)  ~  9(0)] 

V    7TX 


V      7T       (yo  V^  Jo      VI  —  T         J 


From  the  compatibility  conditions  on  the  initial  and  boundary  values,  we  have 


limu(x,0)  =  g(0)  =  limu(0,y)  =  /(0) 

x— >0  y— ♦O 


Therefore 
du 


dy 


(*,y) 


'2jAr 


/-0+ 


7T        [y? 


■7-  f  «r)«-*''»*  -  f    ^Ld 

x/X  Jt=0  Jt  =  0  \/X  —  T 


(35) 


It  is  easy  to  note  from  (35)  that  for  k  =  k0  —  je, 


H^X^ 


=  0. 


y^0+ 
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No  matter  how  small  the  t  is,  we  will  use  the  same  equality  for  the  limiting  case 
of  e  — ►  0.  We  will  evaluate  the  integrals  above  approximately  by  replacing  the 
derivatives  with  differences. 


y 


Ay 


u, 


u. 


tu, 


u 


\ 


1  3 

u     u 


■>    X 


Consider  initial  data  on  the  line  x  =  0.  Let  us  assume  that  this  initial  data  is  known 
on  a  uniform  grid  ym  =  mAy,  m  =  0, 1,  •  •  •.  We  approximate  the  derivative  in  the 
interval  (ym_i,  ym)  by  the  forward  difference  formula 


df{y)      um  -  um_! 


dy 


Ay 


y  e  (ym-uym) 


where  um  =  u(0,ym). 
Then 


Let 


y/x  Jr=o  ^         Ay 

y/X  JT=ym-\ 


—  t  =  \pK  [i  ==>  —j=dr  —  \  -rdfi 
x  \Jx  V  k 
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We  then  have 


y/k/(irx)ym 
y/k/(irx)ym-i 


-l^2l2dil 


m 


"  y/k/{*x)yr> 


•  y/k/(*x)ym-i 


where  F(fi)  =  complex  Fresnel  integral  [8] 


I    \l Vm        -  F\   \     J 

■KX  I  \   V   XX 


(36) 


Consider  now  the  boundary  data  on  the  line  y  =  0.  Assume  that  we  have  a  non- 
uniform grid  0,xi,  x2,  •  •  •  ,xm  —  x.  On  the  interval  (xm_i,£m),  we  approximate  the 
derivative  as 


9t(t)  = 


um  _  um-\ 


xm         xm  —  \ 

where  um  =  u(xm,0).  At  the  origin  we  have  u°  =  u0. 
We  evaluate  the  second  integral  in  (35)  as 

JO      y/X  —  T  m=1   Xm   —  Xm-\   •'im-i    \J X  —  T 


=  £ 


M     um   _  um-\ 


=  \  xm         xm  —  \ 


(-2v/x^7) 


Xm-l 


=    2E 


M    um  -um-1 


771=1 


M 


Xm         xm  —  \ 


[y/X  —  Xm_i  —  y/X  —  Xm  ) 


um  _ um~\ 


-   % ^7^7,^7^ 


(37) 


Using  (36)  and  (37)  in  (35)  we  get 

— -    oo 


Ty^y) 


<2jk 


y=0+ 


71 


'W,/-*-ym  )-F 


m=l 


Ay 


■kxm 


■KXM 


'Um-l 


M 


-2E 


um  _  um-\ 


_j  y/xM  ~  xm  +  y/xM  ~  x 


m-\ 
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Extracting  out  the  m  =  M  term  from  the  latter  and  writing 

duM 


du 
dy 


y-0+ 


we  have 

duM 


dy 


+  M--Z—  =  y/vt 


dy 


Um         ^m  — 1 


Ay 


F 


■KXM 


■Vm]    ~  F 


■kxm 


'Vm-l 


ISjk 


M-l 

E 


um  -  um~l 


*      m  =  \    \/xM  -  Xm  +  \AM  -  lm-1 


(38) 


The  above  equation  is  the  discrete  version  of  a  continuous  boundary  condition  of  the 
form 

-^  +  r(x)u  =  5(x)  (39) 

dy 


where 


and 


r(x)  = 


%k 


■K     yjx  -  XM-] 


(40) 


■to  -  V2;E 

771  =  1 


^771  ^m_] 


Ay 


F,«s«- rF  v^™-1 


M-l 


-Vs?  E 


vT  -  um~l 


7T      f^1>/x^^+>/x-Xm_1     ' 


and 


Fix)  =  re'^^di 

Jo 

is  the  complex  Fresnel  integral  [8]. 


(41) 


(42) 


For  efficient  implementation,  the  PDE  and  the  various  boundary  conditions  will  be 
transformed  into  a  curvilinear  coordinate  system  generated  by  setting  the  lower  ir- 
regular boundary  as  a  constant  curve  curvilinear  coordinate. 


23 


7.      Transformations  to  a  Curvilinear  Coordinate  System 

Consider  the  narrow  angle  parabolic  equation  with  range  dependent  refractive  index 
(aj  ^  0)  given  in  (18) 

1  f  .         d       d2\ 

(2;/:0-ai)  [  dy       dy2 ) 

together  with  the  tropospheric  boundary  condition 

dy 
and  the  impedance  boundary  condition  on  the  irregular  boundary 

uu  4-  C\U  =  0 


<  1-yyv|>*JU~c*       Bawvx^y 


We  will  transform  these  to  a  curvilinear  coordinate  system  (£,77) 


/     ,     / L 


* f  „  c *. i c 


r-r-4  <r-BCL-»; — t     ,    L     ,  1  ,    J  ,  ,   L,   r):Q 


<  <! 


f-  »  -  N 


Physical  Domain 


Computational  Domain 
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We  assume  that  we  have  a  transformation  of  the  following  form 

x  =  x{()  ,  y  =  </(£,//) 

The  various  metrics  needed  in  the  transformed  equations  are  [5] 

9u  =  x\  +  y\  ,  9\2  =  x^Xr,  +  y^yv  . 

We  then  have  with  y/g  —  x^y^  —  xvy^  =  x^yv  ^  0 

ux    =     —{y^ui  -  ytun)  =  {u(yr,  -  u„ft) 

1   /  l  \       Ur> 

uy    =     —-{-xr)ui  +  xiuv)  =  — 

V9  yv 


uyy 


\   d  (uv\       1 


Substituting  into  the  PDE  we  get 


Vr, 


XiVr, 


{u^yr,  -UiVe)  = 


1 


(2jk0  -  d) 


or 


ul 


Letting 


x^a1 


■«  + 


0-2  y-nr) 


Xi 


y^      y?         y; 
xm 


+ 


(2jftb  -  d)         [\yv       y*  /  (2>A;o  -  ai)        yn  J  (2j'A:0  -  ai)yj 


Wr,  + 


jWijU 


&1  = 


(2jAr0  -  ax) 


k  =  a 


a2 


y™ 


6,    = 


S/r,  L\  y?  /  (2j^o-ai) 


+  y* 


(2jko-al)y^ 


we  express  the  PDE  as 


w^  =  6lU  +  62Wr,  +  63u 


*?*] 


(43) 


The  normal  derivative  on  a  77  =  constant  line  is 

/                 .  v              [On              9i2 
uJrj  =  const.)     =     «/  —  un u^ 

V   9  y/9911 


HVri  -  Vixn  yjx\  +  y\(x^yn  -  y^) 


•Mr,    — 


Hxn  +  y^yr 


H 


(44) 


25 


Variation  of  an  arbitrary  vector  r  along  the  77  =  const,  coordinate  is 

H  -  x^x  +  y^y  =  at 
The  unit  vector  along  the  tangent  on  77  =  const,  is  then 


Vi 


Defining 


cos  9    = 


s'm9    = 


>ri 


^l  +  y? 


we  express  (44)  as 


uv  — 


(xv  cos  9  +  y-n  sin  #)t/£ 


y„  cos  0  -  x„  sin  0       ^xf+yfC^  cos  0  -  x„  sin  0) 

With  these  substitutions,  the  boundary  condition  at  the  bottom  boundary  i^  +  cju  = 
0  gets  transformed  to 

(xv  cos  0  +  yn  sin  0) 
it,, ,  it^  +  (y,,  cos  9  —  Xr,  sin  0)citt  =  0     on  r\  =  0  (45) 

v*? +  ^2 

For  x,,  =  0  and  using  Jx\  +  y?  =  x,,/  cos  0,  we  get 


V 

uv sin  9  cos  9u^  +  ciy,,  cos  9u  =  0     @     77  =  0 

x< 


(46) 


Finally  at  the  top  boundary,  we  have 

u„/y„  +  ru  -  s 


or 


u,(^m,  W)  +  r(£m,  W)y„(£m)u(£m,  TV)  =  yv4U,  N) 


(47) 
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8.      Generation  of  Curvilinear  Coordinate  System 


Consider  a  piece-wise  linear  ground  profile  and  a  horizontal  upper  boundary 


t 


-* — . 

— < 

p-^ 1 

A' 

* ' 

i — <= — 
8' 

— « 

^-r1 — *■ 
c 

1 1 

r — <- 

0' 

1 

X 

»          t 

1 

/B\ 

^f^t 

)NX 

k 

c 

Non-Rectangular, 
Uniform  Mesh 

Physical  Domain 


■    '                     i 

/*)  = 

Jructu^ 

A/ 

8' 

c' 

Li-rxgV. ^, 

*  M    — i 

t < 

it 

•--> »— « 

-7 7—* 

6 

Rectangular, 
Uniform  Mesh 

Computational  Domain 


D 


1)z.O 


*JL L L C * 


*z_^ — « — *— « 


^c+t » Yi+i 


9-, 7- 


We  generate  the  (z,y)  coordinates  of  an  interior  point  by  using  a  linear  interpolation. 
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Letting  Ax,  =  (x1+1  -  x,),  Ay,  =  (y,-+1  -  y,),  and  A£  =  (£,+1  -  £,),  we  write 


Ax- 


y=    l- 


A(, 


V 


0  <  77  <  TV 


(48) 


At  any  interior  point,  £,  <  £  <  £,+i>  the  various  metrics  are  evaluated  as 

Ax, 


x(    = 


A6  ' 


xn  =  0 


(49) 


Vi 


"37 


7  A  Ay* 
A6  ' 


?/, 


1 

TV 


2/o  -  Vi  ~  ^T-U  -  CO 


y^f? 


At  the  boundary  points,  we  use  the  central  difference  formulas  to  arrive  at 

xt-i_2  —  xi       Axt-+i  +  Ax, 


*dt=b)    = 


ydt  =  b)  = 


ii+2 -6"  A6+1  +  A6- 


_77_\  Ay,+1  +  Ay, 
"TV;  A{l+i+A6 


(50) 


(51) 


Note  that  the  analytical  expressions  yield  a  discontinuous  value  for  these  derivatives 
at  the  boundary  points.  We  use  the  analytical  expressions  only  to  generate  grid 
points  and  use  the  central  difference  formulas  to  arrive  at  the  derivatives  w.r.t.  £. 
In  other  words,  once  the  grid  points  are  generated  on  the  lines  AA\  BB',  . . . ,  etc., 
we  assume  that  the  space  is  smoothly  connected  through  the  grid  points.  In  the 
numerical  implementation  using  Crank-Nicolson  implicit  scheme  [6],  the  metrics  are 
needed  at  the  midpoint  w.r.t.  £,  i.e.,  at  (  =  £,  +  A£,/2  and  the  interior  point  formulas 
are  applicable.  For  a  uniform  mesh  in  the  computational  domain,  A£,  =  1,  77  =  q, 
9  =  0,1, 2,---,JV 


=    Ax, 


Vi 


1  -  J7.I  ** 


^(9  +  1)  -y&q)  =  - 


Ay, 


Vv     = 


N 


Sfo 


Vi+i  +  Vi 


(52) 
(53) 

(54) 


=       1- 


N 

y.+i  +  y« 


y«+i  +  yi 


+  £y° 


+  qyv  =  yo  +  {q-  N)yr) 


(55) 


so  that  y{q  +  1)  -  y(q)  =  yv. 


28 


9.      Numerical  Implementation  by  Crank-Nicolson  Scheme 

Consider  the  narrow  angle  PE  together  with  the  boundary  condition 

U£  =  biu  +  b2uv  +  b3uvv         (PDE) 

uv  —  I  —  sin  0  cos  0)  Uc  -{■  Ciyv  cos  Ou  =  0     at  tj  =  0         (BC2) 


un  +  ryvu  -  yv  s     at  77  =  N         (BCi) 


with 


uv{0,N)  =  0 
We  would  like  to  implement  the  above  using  a  Crank-Nicolson  implicit  scheme  [6] 


< n * 

i 

n» 1 < 

=-,, ! , 

1 

1 

-HI 1 < 

1 

1 
i 

I 

f>-i 


b-JL 
r      2 


We  use  the  notation  up  =  u(£p,r)q)  =  u(xp,yq) 
The  various  derivatives  assuming  A£  =  A77  =  1  are 

U*(£p-l/2,*7g)      =      <-<_1 

p— 1        p— 1   1     p  p 

MP~  2'V  =  4  ^  +W'+1  ~  (U*~*  +W'"1)] 
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-   n>< 


1 

2  L 


«j+1  -  2U;  +  «*   +  «£}  -  2UJ-1  +  <:J 


Substituting  into  the  PDE,  we  get 


9  ? 


P-l 


Mp~  o'9J    / 
&3  (p-  r,^)     , 

v  22   y  (<+i  -  2<  +  uj.,  +  uj;J  -  2<->  +  <:!) 


Rearranging  the  terms  we  get 


lbi^)u'-.i  +  b3-bi)u'  +  lL-h)uU 


2  V  2 


2  V  2 


6i 


•  63  uj;j  +  1  -  63  +  f  J  tx;-1  +  ^  ( 63  -  f  j  up! 


2 ' 
62 


Now  let 


1/.         1L\H/2 
a    =     2^+262 


,-     (,-| 


8 
p-1/2 


7  =  K63_l 


p-1/2 


We  then  have  for  p  =  1,2,  •  •  •,  q  =  1,2,  •  •  •,  N  -  1 

qU;+1  -  (1  +  /?)<  +  7<-,  =  -«<;]  -  (1  -  /^j-1  -  7<:! 


(56) 


However,  we  will  extend  the  applicability  of  this  equation  over  q  =  0,1,---7V  to 
accomodate  the  derivative  boundary  conditions  on  the  lower  and  upper  boundaries. 
For  q  =  0,  we  have 

au\  -  (1  +  (3)uv0  +  1Up_,  =  -aur1  -  (1  -  0K"1  -  7U-V  (57) 

for  g  =  JV,  we  have 

aupN+1  -  (1  +  /?)<,  +  7"jv-i  =  -«wn~+i  -  (1  ~  /^"n"1  -  7Wn"-i  (58) 
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From  the  BC2  at  the  lower  boundary  we  have 


i — f- 


&=' 


? t« — X — » 7 7-      |:0 


+ 


V- 


--I 


U-l  u\  —  W_j 


+ 


\P~l/2 

—  sin  ^  cos  ^  J  \Uq  —  Uq~) 


+  (c!yv  cos  0)1 


p-l/2  /  uo  ~  uo 


=  0 


p-1/2 


Multiplying  this  with  4(7)0         an<^  adding  to  (57) 


(a  +  7)t^-  (l+/?-27 
=  -(a  +  7)"i_1-  I  1-0  +  27 


,S/r 


Cj  j/„  cos  0  —  2 —  sin  0  cos# 


p 


V 

Ciyv  cos  9  +  2—  sin  6  cos  0 


p-i 


Letting 


a 


a  +  7 


/?'    =    /5  -  27^  cos  0  ( Ci  - 
0"    -    £-27yr;cos0fCl  + 


2sin0' 
**     . 

2sinfl' 
xi    , 


we  may  write  the  above  equation  as 

a'ul  -  (1  +  /?X  =  -a'u\-1  -  (1  -  /?")^_1 


From  the  BCi  at  the  top  boundary  we  have 


J. L. 


I>- 


X    '  * 


—  Nl+I 

-*    N 

-  N- 


(59) 
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p-1 

uyv+i  ~ 


p-i  j. 

N-l     ,    uN+l  - 


V 
N-l 


+ 


rv  +  rp-i      up   +  u 


-^ 


p-1 

AT  T"  U/V 


SV 


SP  +  S 


p-1 


2  2 

Multiplying  this  with  — 4(a)^"       and  adding  to  (58),  we  get 

-  (l  +  /?  +  [rp  +  *•"-%„<»)  txj,  +  (7  +  4a)upv_1 

-  -  (l  -  0  -  [rp  +  r""1]^)  u^1  -  (7  +  4a)WAT_11  -  2yr,(5p  +  s^1) 
Letting 

A    =    p  +  y^r"  +  rp~'} 
7'     =    7  +  4q 

we  write  the  above  equation  as 

(1  +  A)<  -  iu*N_,  =  (1  -  AX"1  +  7'<-_\  +  2j,„(r>  +  r?"1) 


(60) 


We  augment  the  equations  given  in  (56)  for  q  =  1 , 2, . . . ,  N  —  1  with  (59)  and  (60)  for 
q  =  0  and  N,  respectively  to  define  for  q  =  0, 1, 2, . . . ,  TV.  The  system  of  equations  so 
defined  can  be  expressed  as  a  matrix  equation  of  the  form 


X    X 
XXX 
XXX 


X    X 


lN 


X    X 
XXX 
XXX 


X    X 


V 


u 


p-1  "I 

0 

p-1 


.p-1 

'.V 


(61) 


where  X  denotes  a  non-zero  entry.  The  tridiagonal  matrix  on  the  left  hand  side  of 
(61)  can  be  inverted  efficiently  to  yield  a  solution  on  line  £  =  £p  in  terms  of  the  field 
values  on  the  line  £  =  £P-i-  Equation  (61)  can  be  used  to  march  forward  in  range 
starting  from  initial  data  specified  on  £  =  0. 
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